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Abstract 

For the 2D Euler equations and related models of geophysical flows, 
minima of energy-Casimir variational problems are stable steady states 
of the equations (Arnol'd theorems). The same variational problems 
also describe sets of statistical equilibria of the equations. This paper 
uses Lyapunov-Schmidt reduction in order to study the bifurcation di- 
agrams for these variational problems, in the limit of small energy or, 
equivalently, of small departure from quadratic Casimir functionals. 
We show a generic occurrence of phase transitions, either continuous 
or discontinuous. We derive the type of phase transitions for any do- 
main geometry and any model analogous to the 2D Euler equations. 
The bifurcations depend crucially on 04, the quartic coefficient in the 
Taylor expansion of the Casimir functional around its minima. Note 
that ci4 can be related to the fourth moment of the vorticity in the 
statistical mechanics framework. A tricritical point (bifurcation from 
a continuous to a discontinuous phase transition) often occurs when 
04 changes sign. The bifurcations depend also on possible constraints 
on the variational problems (circulation, energy). These results show 
that the analytical results obtained with quadratic Casimir functionals 
by several authors are non-generic (not robust to a small change in the 
parameters) . 



1 Introduction 

Flows which are turbulent and two-dimensional are remarkable for two rea- 
sons: First, they self-organize into large-scale coherent structures and, sec- 
ond, they often display a bistable behaviour. Such large-scale structures 
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(monopoles, dipoles, parallel flows) are analogous to geophysical cyclones, 
anticyclones, and jets in the oceans and atmospheres [5]. The motion of atmo- 
spheres and oceans is almost two-dimensional indeed, because of the following 
three characteristics: i) the fluid layer has a small vertical-to-horizontal as- 
pect ratio, ii) the fluid layer is subject to the Coriolis force, which dominates 
the viscous frictional forces, and iii) the water column is stably stratified over 
large scales, constraining motions to be horizontal. In the laboratory, these 
characteristics can be obtained by rotating a (shallow) cylinder filled up with 
water and using a forcing mechanism. Typically, observations made on these 
experimental systems can shed light on atmospheric and oceanic phenomena. 

The natural equations governing this type of motion are the Navier- 
Stokes equations in two dimensions. It should be noted that the large scales 
of geophysical flows are highly turbulent. Indeed, scale analysis shows that 
the motion of the large scales is dominated by the advective (also called 
inertial) term; forcing and dissipation terms are small with respect to the 
inertial term. We say that the flows self-organize, precisely because the large- 
scale structures are not at all determined (say, linearly) by some external 
forcing. This self-organization of the large scales is specific to 2D turbulence 
[10]. Unlike in 3D turbulence, there is no direct energy cascade (towards 
small scales) but there are an inverse energy cascade (towards largest scales) 
and a direct enstrophy cascade. This difference stems from the conversation 
of different quantities for the conservative dynamics — then, the 2D Navier- 
Stokes equations reduce to the 2D Euler equations. If topography is included 
in the model, we have instead the (inviscid) barotropic quasi-geostrophic 
equations. In this inertial limit, the attractors of the dynamics are expected 
to be found near a set of steady states of the inviscid equations. 

We have the transport of a scalar quantity q by an incompressible two- 
dimensional velocity; for the 2D Euler equations, q is the vorticity and, for the 
barotropic quasi-geostrophic equations, q is the potential vorticity. We wish 
to predict the final state (s) of the system. Any state verifying a functional re- 
lationship between (potential) vorticity and streamfunction is a steady state. 
Thus, there are an infinity of steady states. How can we determine which ones 
are stable? Since q is a field, the system has an infinite number of degrees of 
freedom (continuous system). A deterministic approach would be unrealistic. 
Then, we turn to statistical mechanics. Rather than describing fine-grained 
structures (exact fields), equilibrium statistical theories of two-dimensional 
turbulent flows predict — assuming ergodicity — the final organization of the 
flow at a coarse-grained level: a mixing entropy is maximized under the con- 
straints that all the flow invariants be conserved [12, 13, 14]. There are an 
infinity of invariants, namely, the energy and the Casimirs; a Casimir is any 
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functional of the (potential) vorticity. 

The Miller-Robert-Sommeria theory (MRS theory, for short) predicts 
statistical equilibria in terms of a functional relationship between (poten- 
tial) vorticity and streamfunction. We want to determine the large scales 
of (quasi) two-dimensional turbulence as equilibria of the inviscid equations. 
The analytical computation of MRS equilibrium states would be a difficult 
task though: it would be about solving a variational problem involving an 
infinite number of constraints. In this paper, we present analytical and nu- 
merical computations of phase diagrams for a large class of equilibrium states, 
obtained from simpler variational problems. Phase transitions are very im- 
portant to study, since they are associated with major physical changes (in 
large-scale flow structures, as far as we are concerned) in the system under 
consideration. For instance, flows will have their structure change as they un- 
dergo phase transitions. It is important to know whether these are first-order 
(discontinuous) or second-order (continuous). 

Indeed, simpler variational problems (taking into account only a few con- 
straints) were shown to give access to some classes of MRS equilibria [2] . For 
instance, one such class is the one for which q = f{ip) is linear (or affine). 
An example of using statistical mechanics for predicting and describing real 
turbulent flows can be found in [4] and references therein. Bifurcations be- 
tween stable steady solutions of 2D Euler are found to occur when varying 
the domain shape, the nonlinearity of f(ip), or the energy. This suggests that 
a general theory of phase transitions for 2D and geophysical flows should be 
looked for — it is not available at the present day. Only instances of such 
phase transitions have been reported in the literature. Note that key results 
regarding statistical ensemble inequivalence, encompassing the case of a non- 
linear equation q = f(ip), were presented in [9]. In this paper, we present new 
analytical results on phase transitions related to the nonlinearity oi f(ijj) . We 
obtain a complete theory of phase diagrams for two-dimensional turbulence 
equilibria and steady states in the low-energy limit. 

The simpler variational problem we consider writes 

C S (E, r) = min { jf s(q) | £[q] = E , r[q] = T j . (1) 

The function s(q) is assumed strictly convex. In thermodynamics, the micro- 
canonical problem is a two-constraint variational problem where the thermo- 
dynamical potential to be maximized is called the entropy. We can draw an 
analogy with (1), where our Casimir functional J v s(q) acts as the opposite 
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of an entropy. We give the expressions 1 of 

1 f 

£[q] = — / ijj{q — h) , the kinetic energy, and 

2 Jv 

r[o\ = / Q i the circulation. 
Jv 

So we call (1) microcanonical, in analogy with usual thermodynamics. Note 
that this variational problem corresponds to (CVP) in [2] (see this reference 
about the relationship between the solutions to our variational problem and 
the actual MRS statistical equilibria). For given values of the constraints E 
and T, the q fields solving (1) are microcanonically stable equilibria. This is 
a sufficient condition for their dynamical stability [1]. Indeed, let us consider 
a functional which is conserved by the dynamics. This functional can be a 
linear combination of a Casimir and of the energy ('energy-Casimir func- 
tional'). The point is the following: if the system lies at a nondegenerate 
extremum of this invariant, then it cannot go away from this point. 

The paper is organized as follows: In section 2, we define the quantities 
and notions in use throughout the work and we give the general results. It 
appears that phase transitions can be characterized through the bifurcation 
analysis of scalar equations, the latter acting as normal forms. The technical 
derivation of the results is given in the various appendices. In section 3, 
we apply the general results to a particular case, in a rectangular domain. 
Equilibria are computed numerically using the pseudo-arclength continuation 
method. Finally, in section 4, we suggest some physical applications, offering 
a clear motivation for this theoretical work. 



2 Definitions and general results 

The system we consider is that of the barotropic quasi-geostrophic equations, 
which model the 2D dynamics of one oceanic or atmospheric layer: 

d t q + u-Vq = 0; u = e z x Vip ; q = Atfj + h (2) 

where u denotes the (two-dimensional) velocity field, ip the streamfunction 
(defined up to a constant), q the potential vorticity (in vorticity units), and h 
an equivalent topography. The boundary condition is ip = on &D, where D 
is a simply connected domain in two dimensions. The natural scalar product 
for the fields at play is denoted by (qiq?) := J v qiq2- 

1 See (2) and section 2 for a definition of the fields ip and h. 
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The inviscid dynamics (2) corresponds to a limit of infinite Reynolds num- 
ber. Although the number of degrees of freedom is infinite in a turbulent flow, 
the formation of large-scale structures indicates that just a few effective de- 
grees of freedom should be enough to characterize the flow. In this paper, we 
describe a class of steady states of (2), and the phase transitions which they 
undergo, through scalar bifurcation equations. The stability of these steady 
states can be established statistically (thermodynamically), implying dynam- 
ical stability (Arnol'd's theorems). In the equilibrium statistical-mechanical 
context, we deal with phase transitions. In the dynamical system context, 
we deal with bifurcations. Here, we use technical tools of applied bifurca- 
tion theory, namely, Lyapunov-Schmidt reduction, to characterize the phase 
transitions. Hence, we can determine the continuous or discontinuous nature 
of phase transitions in a general framework. 



2.1 Relaxed variational problems 

To compute statistical equilibria, which are, again, particular steady states 
of the dynamics (2), we solve the microcanonical variational problem (1), as 
announced in the introduction. In this paper, we restrict our attention to 
even functions s(q). Indeed, there are many situations where the q i— > —q 
symmetry applies. If q is a solution to (2), then — q is also a solution to (2). 
In real flows, the q y —q symmetry could be broken by a nonsymmetric 
forcing or by a nonsymmetric initial distribution of (potential) vorticity. Say 
that s can be written as the expansion 

n>2 

Assuming that the Lagrange multiplier rule applies (q regular enough), there 
exists a couple (/3,7) € M 2 such that solutions of (1) are stationary points of 

g[q]= I s (q) + /3E[q}+jr[q}. (4) 
Jv 

We call this functional the Gibbs free energy, in analogy with usual thermo- 
dynamics. The variational problem dual to (1), i.e., 



G(/3,7) = min<j£[g] = J s(q) + f3£[q] + jr[q] j , (5) 

is referred to as the grand canonical variational problem. Because it is relaxed 
(unconstrained), it is more easily tractable. 
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We shall consider another energy-Casimir variational problem, namely, 
the canonical variational problem: 

F{p, Y) = min ^F[q] = jf s(q) + 0S[q] \ r[q] = T J . (6) 

It is the problem of minimizing the Helmholtz free energy with fixed circu- 
lation r. 

2.2 Ensemble inequivalence 

For given values of the constraints E and T, the q fields solving (1) are 
statistical equilibria. As introduced in the previous subsection, (3 and 7 
are the Lagrange multipliers associated with the energy and circulation con- 
straints, respectively. For all couples (/3,7), minima G(/3,7) are also minima 
C s (E(/3, 7), r(/3, 7)). But some minima C S (E,T) may correspond to station- 
ary points of (4) which are not minima of (4). These are classical results 
(see any textbook on convex optimization). When E(^,j) and r(/3, 7) do 
not span their entire accessible range (E € M + , r G 1) as (/?, 7) is varied, the 
microcanonical ensemble and the (dual) grand canonical ensemble are said 
to be inequivalent. Then, some microcanonical solutions are not obtained as 
grand canonical solutions. 

This feature is typical of long-range interacting systems. We use the 
term 'long-range interactions' as in, for instance, [3] and [6]: for a system 
in space dimension D, the interaction potential between particles separated 
by a distance r goes like r~ a , as r — > 00, with a < D. The interaction is 
'non-integrable'. From the expression of kinetic energy for 2D Euler, the 
coupling between vorticity at point r and vorticity at point r' appears to 
be logarithmic, hence not integrable. Thus, the vorticity at a given point 
is coupled with the vorticity of any other point of the domain, not only of 
neighbouring points. In addition to 2D turbulence, long-range interacting 
systems include self-gravitating systems in astrophysics and some models in 
plasma physics. In short-range interacting systems, the different statistical 
ensembles are used interchangeably, because they are usually equivalent. 

Let us illustrate the idea of ensemble inequivalence with a schematic pic- 
ture. For the sake of simplicity, let us discard the circulation constraint. 
The microcanonical solutions are described by C S (E). If the caloric curve 
/3(E) = —C' S (E) is monotonically decreasing, i.e., if C S (E) is convex, the 
microcanonical solutions can all be obtained as canonical solutions: the two 
ensembles are equivalent. If the caloric curve is increasing over a certain 
range (negative specific heat, in thermodynamics terms), there is a range of 
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ensemble inequivalence. The canonical ensemble is equivalent to the micro- 
canonical one only over the range for which C S (E) coincides with its convex 
envelope (range E > E c on Figure 1): solutions C S (E) are solutions F(f3). 
C S (E) has an inflexion point at E = E C2 (we shall use this notation in sub- 
section 2.4). It is a canonical spinodal point. We refer the reader to [3] for a 
systematic classification of all these singularities. 




Figure 1: C S {E) showing a range of 'canonical stability' for E > E c : minima 
C S (E) are minima F((3); a range of 'canonical metastability' for E € [E C2 ,E C ]: 
minima C S {E) can be obtained as local minima of the canonical functional; a 
range of 'canonical instability' for E € ]0, E C2 [: minima C S {E) can be obtained as 
local maxima of the canonical functional. 



2.3 Poincare inequality 

The Poincare inequality comes in handy to establish a sufficient condition 
for convexity; it is natural to begin with the study of the convexity of Q[q] 
(4). Indeed, it is readily noted that if Q[q] is strictly convex, it has a unique 
stationary point, which is then the (unique) solution of (1). Since r[q] is a 
linear form, it is sufficient to investigate the convexity of the Helmholtz free 
energy functional F[q] = f v s(q) +j3S[q], Since S[q] is convex, F[q] is strictly 
convex if > 0. 
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If /3 < 0, we need to study the sign of the second-order variation of J 7 , 
denoted by b 2 T , and defined through F\q + Sq] — JF[q\ = 8JF[q\ + \5 2 J r [q] + 
o(5q 2 ). We get 

5 2 T[q] = [ s"(q)5q 2 - (5 [ 5^5q. (7) 
Jv Jv 

We make use of the Poincare inequality: 

f > ~ I Sq 2 , 

Jv M Jv 

where — Ai < is the greatest (smallest, in absolute value) eigenvalue of the 
Laplacian on T>. Indeed (let us recall the classical proof), introducing the 
orthonormal Laplacian eigenbasis {ej(r)}j>i, i.e., 

Aej(r) = -Aiei(r), / = < Ai < A 2 < . . . , (8) 

Jv 

with Ci = on dT> for alH > 1. All fields may be decomposed in this basis: 

i 

5q(r) = 5(A^(r) + h(r)) = 5A^(r) = -£ Wi^(r) 

= ^% ei (r). (9) 

i 

Therefore, 
So 

5 2 J-[q] > jf (V(<?) + V > + jf V, (10) 

for (3 < 0, where s'^ := min re 2){min 9 s"(q(r))}. If (3 > —s'^Xi, J 7 is strictly 
convex, and so is Q. There is a unique solution to (5) and, hence, a unique 
solution to (1). 

The conditions ft > and —s'^Xi < < are the hypotheses for the first 
and second Arnol'd theorems, respectively, on Lyapunov stability. In both 
cases, the sufficient condition is that 5 2 J r be positive-definite [11]. We can 
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conclude that for (3 > — s'' t Ai, the grand canonical ensemble is equivalent 
to the microcanonical ensemble. In the grand canonical ensemble, phase 
transitions may occur only for (3 < —s'^Xi, where solutions to (5) may cease 
to be unique or cease to exist. 

The stationary points of Q (4) are the q fields for which the first-order 
variation of Q vanish, i.e., 

s'(g)-/3V + 7 = 0. (11) 

Since s(q) is strictly convex, s'(q) is strictly increasing, so its inverse (s')~ 1 (q) 
is well-defined (and strictly increasing). We have 

q=(sT 1 (^-l)- 

From (3), the Taylor expansion of (s') _1 around reads (s')~ 1 (x) = x + 
CL4X 3 + o(x 4 ). Then, the term in CI4 is the lowest-order nonlinear contribution 
to (s') _1 (x). 

2.4 Phase diagram for 7 = 

It is more straightforward to study a symmetric problem first and, afterwards, 
to study the effect of breaking the symmetry. Therefore, we begin with 
the case 7 = so that (5) is symmetric with respect to g 4 —q. The 
corresponding constrained variational problem is the grand microcanonical 
variational problem with 7 = 0, i.e., the minimization of J v s(q) with fixed 
energy. We find that the grand canonical ensemble with 7 = is equivalent 
to the grand microcanonical (only energy-constrained) ensemble if 04 < 0. It 
is not the case if 04 > 0. 

We have denoted the first (largest-scale) Laplacian eigenmode by e\. As 
long as the topography field h is orthogonal to e±, we find the following 
results, for the grand canonical ensemble with 7 = 0: 

• for 04 < 0, there is a second-order phase transition at = —Xi. the 
solution goes continuously from a trivial state (zero energy, uniform 
vorticity) to a state dominated by e\\ 

• for 04 > 0, Gt4 small enough, there is a first-order phase transition at 
P = PciflA) £ ] — -^l) ~^i s m[' the solution goes discontinuously from a 
trivial state (E = 0) to a state dominated by e\ (E = E c {ai) > 0). 
The energy range accessible by grand canonical solutions (with 7 = 0) 
displays a gap ]0, E c (a^)[. 
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Systems with symmetry display a richer phenomenology of phase transitions, 
especially regarding second-order phase transitions [3]. So it is not surprising 
to find a second-order phase transition line here. 

In the grand microcanonical ensemble with 7 = 0, we find that 

• there is no phase transition at low energy (we cannot tell what happens 
at high energy); 

• at nonzero low energy, the solution is a state dominated by e\\ 

• for a 4 > 0, states of lowest energy (E G [0, E C2 {a^]) have negative 
specific heat. 

What is the method for deriving these results? We explain it qualitatively 
here and give the technical details in appendix 5. When solving (1), the 
quadratic part of s comes into play at lowest (linear) order in E [4] . Therefore, 
in the low-energy limit, it is the dominant contribution. Also, at lowest order, 
the solution is along e±, the largest-scale eigenmode. The next order brings 
into play the small parameter 04, referred to as the nonlinearity, for short. 
We may always write q = Aei + q' with A G R and q' orthogonal to ex- We see 
q' as a perturbation to the lowest-order solution ±Aex and assume it admits 
an asymptotic expansion in (powers of) A. This will lead to an asymptotic 
expansion in A for the Gibbs free energy, i.e., a normal form describing the 
phase transitions in a neighbourhood of a 4 = 0. The idea is to minimize Q 
with respect to q' first, then with respect to A. We expect the symmetries 
at play to show in this normal form. 

Thus, we have reduced the infinite-dimensional variational problem — or 
equation for the stationary points (11) — to a scalar equation, the bifurcation 
equation [8]. This is called Lyapunov-Schmidt reduction. We have solved 
the unconstrained (or grand canonical) variational problem in a symmetric 
case. We find a tricritical point. 

2.5 Tricritical point 

A tricritical point is a point where a second-order phase transition meets 
a first-order one. We have predicted the phase diagrams in the vicinity of 
(/?, 04) = (— Ai, 0). For 7 = and hi = 0, the normal form (17) is 

Go(A) = \ (l + £) A*-^(et)A* + (A% 

We can readily relate this expression to the normal form s aj b(m) = — m 6 — 
36m 4 /2 — 3am 2 . This normal form is used in the context of constrained 
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variational problems in [3]. Note that s a ^(m) is to be maximized, and hence, 
solutions are maximizers there. Then, the identification of coefficients is to be 
done between s a ^ and —Go- The typical behavior of s a ,b and the associated 
transition lines are shown on Fig. 6 of this reference, reproduced below (our 
figure 2). 



b 




Figure 2: The 'canonical' tricritical point is at (a, b) = (0,0). The curve (4a = 
b 2 ,b < 0) corresponds to the appearance of three local maxima. The bold curve 
(16a = 3b 2 , b < 0) is a first-order phase transition line. The bold-dashed curve is 
a second-order phase transition line. Here, 'canonical' simply refers to a relaxed 
ensemble with respect to a constrained one. Figure from [3]. 

If a > and b > 0, then s a ^ is concave and there is only one maximizer, 
namely m — 0. We can see that m = is always a critical point. The other 
possible stationary points are such that m 4 + bm 2 + a = . For b > , a 
pair of maxima appears as a becomes negative, hence the second-order phase 
transition at (a = 0, b > 0). A pair of minima and a pair of (local) maxima 
appear as \b\ > 2y/a in the quarter plane (a > 0,b < 0). There is a first- 
order phase transition when these maxima reach s a ^{m = 0) = as a and b 
decrease. It is found to occur for 16a = 3b 2 . 

Parameters a and b are identified with (1 + ^-)/6 and — a 4 (e 4 )/6 respec- 
tively. Therefore, (0,0,4) = (— Aj.,0) is a tricritical point in the grand canon- 
ical ensemble with 7 = 0. The normal form sketched in the various areas 
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of the phase diagram (Figure 2) should be identified with the opposite of 
Go(A). Stationary points near A = are found to be A = for all (/3, 04); 
in addition, 

are stationary points when (3 + Ai and a 4 have the same sign. 

- For > — Ai, Aq is a local minimum. 

- For /3 < — Ai and 04 < 0, {A±\ are local minima originating from 
symmetry-breaking: there is a second-order phase transition at (/3 = 
— Ai, 04 < 0). 

- For /3 > — Ai and 04 > 0, {A±} are local maxima. Then, minima far 
away from A = have to exist, for Go (A) has a lower bound, owing to 
the convexity of s(q). These, say, 'nonlocal' minima cannot be obtained 
perturbatively. 

- For (3 < — Ai and 04 > 0, A is a local maxima; it is the only sta- 
tionary point obtained perturbatively. Solutions have to be the above- 
mentioned nonlocal minima. So there has to be a first-order phase 
transition at > — Ai, where the solution jumps from A$ to the 'non- 
local' minima. 

Since we also know (see subsection 2.3) that Aq is the only solution for 
(3 > — Ais^, the first-order phase transition is a line /3 c (a 4 > 0) such that 
/3 c (a 4 ) e] - Ai, -Ais^(a 4 )[. 

2.6 Phase diagram for constant circulation 

Now we solve the canonical variational problem: the energy constraint is 
relaxed, the circulation is fixed at a low value. We find interesting phase 
transitions, where the flow structure completely changes. For elongated rect- 
angular domains (aspect ratio r > r c ), we recover the showing up of a dipolar 
structure (contribution from mode e[), while for square-like domains (r < r c ), 
we recover that of a central monopole with counter-circulating cells at the 
corners (contribution from mode e*) [7, 16]. The novelty here is to distin- 
guish between a first-order transition and a second-order one, depending on 
the sign of the nonlinearity in q = f(ip), at zero circulation. First of all, we 
restrict our study to the case of zero circulation (r = 0), bringing symmetry 
to our system. As noted earlier, systems with symmetries are well known 
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to display a richer phenomenology of phase transitions. We obtain phase 
diagrams with tricritical points, again. Thus, results in the microcanonical 
ensemble (with T = 0) can be deduced the same way as in subsection 2.4: 
for 04 < 0, no singularity of C S (E); for > 0, canonical spinodal point, 
negative specific heat for the lowest-energy states. 

In the linear case (04 = 0), the constrained canonical problem was trans- 
formed into a tractable equivalent unconstrained problem [16]. We use the 
same trick, as detailed in appendix 6. In appendix 8, we detail the com- 
putation of the T = solutions. Because we linearize (11), it is natural to 
recover the same critical values (collectively denoted by A c ) and neutral di- 
rections (collectively denoted by e c ) as in the linear case (a^ = 0), which was 
investigated by [7, 16]. 

At small but nonzero circulation, we lose the second-order phase tran- 
sition to symmetry-breaking, but then we have metastable states (of which 
stability can be made as close as wanted to that of the equilibrium, as the 
circulation tends to zero). In the square-like case, we can be in the presence 
of three qualitatively different states (stable or metastable). 

Let us consider the phase space 04). Right to the first-order phase 
transition line, the solution is a weak monopole (the amplitude A of e c is very 
close to 0). As the first-order phase transition line is crossed, \A\ jumps to a 
larger value, giving a different structure to the solution flow. For example, in 
case ii), the transition to a dipolar contribution is abrupt in the upper half- 
plane, while it is smooth in the lower half-plane, with a canonical metastable 
state showing up (local minimum). The reader is referred to Figure 11. 

Figure 3 shows a schematic phase diagram for case i). Equilibrium states 
of the left-hand side have different topologies, depending on the relative con- 
tributions of the monopole and e*. The contribution of the monopole is 
determined by |T|, that of e* by For certain values of T, there is a re- 
gion in the left-hand-side neighborhood of (/3,a4) = (— A*,0) where the two 
contributions have the same order of magnitude, yielding a tripolar structure 
for the equilibrium states. At large \A\ (i.e., very negative /3, at given 04), 
only e* contributes to the structure of the equilibrium states. 

We wish to emphasize that the canonical ensemble may be relevant to geo- 
physical applications, since the two regimes of known bistable systems have 
different energies. The area of phase diagram near the discontinuous tran- 
sition should be that of interest, when investigating stochastically induced 
transitions. 
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Figure 3: Schematic phase diagram in the canonical ensemble for a rectangle 
of aspect ratio 1.1 (case i)). The blue curve represents the first-order phase 
transition line. This first-order phase transition line ends at a (microcanon- 
ical) critical point, which is located above and close to the point (— A*,0). 
Insets show vorticity fields. The color scale shows negative (resp. positive) 
values in blue (resp. red); the black contours are ten iso- vorticity lines on 
each plot. From left to right: equilibrium dominated by e*; equilibrium con- 
sisting of equivalent contributions from e* and the low-circulation monopole; 
low-circulation monopole. 
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3 Example: Rectangular domain 



In this section, we apply our general results to the (simple) case of a rectan- 
gular domain of area unity: T> = {(x,y) E [0, r 1 / 2 ] x [0,r -1 / 2 ]} with r > 1. 
We choose a function s(q) such that 



with ci4 G [—1/3, 1/6] so that s(q) is convex, as required. Bound = —1/3 
corresponds to q = t&nh.(/3if}) (two-level vorticity distribution {±1} in the 
MRS theory), while bound = 1/6 corresponds to q = sinh(/3^) (three-level 
vorticity distribution {±1, 0} in the MRS theory). We have = 04/4 — 7/60. 
We take h = 0. 

The corresponding steady states are computated numerically by a method 
of continuation, namely, pseudo-arclength continuation. Pseudo-arclength 
continuation is well-suited for computing solution branches which undergo 
bifurcations. Our continuation parameters are the control parameters, /3 and 
CL4. 

3.1 Solutions for 7 = 

We begin with 7 = and r = 1 (square domain). We solve Aip = (s / ) _1 (/3'0) 
(11) in if;, that is, we compute the stationary points of Q§. Thanks to the 
parity symmetry, we may restrict our study to the domain A > 0. For a 
given 04 > 0, A + (12) is the local maximum. If we increase /3 from — up 
to /3 C2 , we can bifurcate into the 'nonlocal' minimum of Gq(A), as represented 
on Figure 4. Thinking of A as an order parameter, there is a fold bifurcation 
at (5 = (3 C2 (a A ). 

In a square domain V, A x = 27r 2 w 19.7392 and ei(x, y) = 2 sin(7nr) sin(7n/). 
We start at (/3, 04) = (— Ai + 0.006, 0.015) with solution guess 



Let us denote by A comp the scalar product of the (computed) solution q with 
mode e±. \A comp — A + \ must scale like A+. We check that we caught the 
proper solution branch by verifying this scaling relation. Figure 5 shows 
^4comp as a function of /3, displaying the expected fold bifurcation. 

We show the value of Go as a function of /3 on Figure 6. The first-order 
phase transition (/3 = /3 C ) is found as Go {A ^ 0) vanishes. We compute the 
line j3 c {ai > 0) using continuation on (3 and on 04. Just like 16a = 36 2 is the 
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Figure 4: Schematic representation of Go (A) for a given value of > and 
different values of f3, in increasing order from (3 < —X\ (left) to (3 C2 (right). The 
red bullets (stationary points) are numbered according to the path taken by the 
continuation computation: '1' and '2' are A + ] at '3' we bifurcate into the 'nonlocal' 
minimum of Gq(A). 




-19.735 -19.734 -19.733 -19.732 -19.731 -19.73 -19.729 





Figure 5: A{fi) for a 4 = 0.015. 
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-19.74 -19.73 -19.72 -19.71 



Figure 6: Gq(/3) for = 0.030. Local maxima are in blue; local minima in 
magenta; global minima in red. 
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first-order phase transition for the normal form s a> 6 (Figure 2), we recover 
the scaling 

+ AI. 

on the first-order phase transition line, as shown Figure 7. 

The phase diagram in the grand canonical ensemble (7 = 0) is shown 
Figure 8. Figure 9, we show the caloric curve for a positive value of a^. 



3.2 Triple point for 7 ^ 

We show the effect of having topography not orthogonal to the largest-scale 
mode e\ and of having 7 7^ 0: this is the general case for the grand canonical 
problem. Then, the normal form reads 

G{A) = (< ei >7 - A + \ (l + £) A- - |<e^ 4 + 

+ 0{A 3 1 ,A 6 ,A 3 1 3 ,A 4 1 2 ). 

We see that the effect is that of breaking the A h-> —A symmetry of G(A), 
which is a normal form for the grand canonical potential (to be minimized). 
We may take hi = without loss of generality, because the effect of hi ^ 
is qualitatively encompassed by 7 ^ 0. 

Since the second-order phase transition we had originated from the A y 
—A symmetry, we lose it in the general case 7 ^ 0. Therefore, the tricritical 
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10 3 10"' 

l 



10" 



Figure 7: log(a4) as a function of log(l + /3 c («4)Mi)- The cyan and blue points, 
almost superimposing, are upper and lower bounds respectively, for the first-order 
phase transition line. The sets of points are well fitted by a straight line of slope 
1/2. 

point is lost. We are left with a critical point, when the first-order phase 
transition survives. It does so for small enough |7|. It simply gets shifted in 
phase diagram 04): now, (3 C depends on both CI4 and 7. We illustrate this, 
at given small a 4 > 0, in figure 8. In the grand canonical ensemble, we have 
a triple point in phase diagram (/?, 7). 



3.3 Solutions for T 7^ and r = 2 

Figure 11, we show the computed solution for low but nonzero circulation 
(canonical ensemble) and r = 2 — elongated rectangle, case ii). 
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Figure 8: Phase diagram in the vicinity of the (grand canonical, 7 = 0) tri- 
critical point — black dot at (— Ai, 0), where the second-order phase transition 
— green dashed line at (— Ai,a4 < 0) — and the first-order phase transition 
— line (/3 c (a4),a 4 > 0) between the light blue dot series and the dark blue 
dot series — meet. Insets show vorticity fields at (f3 c — (0.006 ± 0.001), 0.030) 
and at (— Ai — (0.006 ±0.001), —0.030); color scale ranges from to 0.6 (from 
blue to red). 



-19.71 



&q -19.72 

CO 



CO 



-19.73 



-Ai 




0.5 



1.5 



2.5 



Figure 9: Caloric curve for 04 = 0.030. 
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. x 10"' 



-1.5 1 

-19.738 



A< 



o-s - A > 



Pc Pc 



Figure 10: Phase diagram in dual space (/3,j) at a 4 = 0.015, display- 
ing a first-order transition line at 7 = up to /3 C > — Ai, splitting into 
two first-order transition lines (7 y —7 symmetry) for /3 G ]/3 c ,/3 C2 [ (in- 
sets below the phase diagram show sketches of a sixth-order normal form 
for G(A), when > 0: different curves on each diagram correspond to 
7 = {-0.001, -0.01, -0.02, -0.03, (-0.05)} from top to bottom, when looked 
at in domain A > ; lhs is for — Ai < < (3 C , rhs is for (3 C < f3 < /3 C2 ). Values 
of A for the solution states are shown in the different regions of dual space. 
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Figure 11: Phase diagram in the canonical ensemble (circulation equal to 
0.01) for a rectangle of aspect ratio 2 (case ii)). The blue curve plots = 
0.13^/Ai + 0; the 0.13 prefactor was chosen so as to fit the three first-order 
phase transition points computed from numerical continuation. This first- 
order phase transition line ends at a second-order phase transition point 
— green dot at (— A' l7 0). Insets show vorticity fields at (/3~, 0.030) and at 
—0.030); color scale ranges from —0.5 to 0.5 (from blue to red); the 
black contours are ten iso- vorticity lines on each plot. 
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4 Physical applications 



First, let us mention that there is a theoretical interest for a classification of 
phase transitions. Notably, two-dimensional flows are long-range interacting 
systems (see subsection 2.2). The nice thing about these systems is that 
theoretical results for one of them is relevant and useful to the others. So 
our general results can extend to other systems with long-range interactions. 

Let us come back to the geophysical motivation though. The quasi- 
geostrophic equations serve as a simple model for the motion of atmospheric 
or oceanic flows. The inviscid quasi-geostrophic model is thus an appropri- 
ate model for geophysical and experimental flows on time scales much less 
than the dissipation time scale, but large enough for turbulent mixing to 
have operated as much as allowed by the constraint of energy conservation. 
For instance, such an equilibrium approach led to a successful description of 
the self-organization of Jupiter's atmosphere; especially, it led to a realistic 
model for Jupiter's Great Red Spot and other vortices [5]. 

Recently, it was shown that ocean currents, such as the Kuroshio (in 
the north Pacific ocean, off Japan) or the Gulf Stream, may also be under- 
stood as equilibria of the inertial dynamics, in very simple ocean models [17]. 
Naturally, this conservative theory ignores all effects due to forcing and dis- 
sipation, which are present in any real flow. Still, a recent work showed that 
the inertial description of equilibria is fundamental and relevant even in the 
presence of forcing and dissipation [4]. Fluctuations are responsible for the 
phenomenon of transitions between two equilibria (bistability). 

The barotropic quasi-geostrophic model is also relevant to the description 
of experimental flows, such as the approximation of fluid dynamics when 
three-dimensional motion is constrained by a strong transverse field (e.g., 
rotation) or takes place in geometries of small (vertical-to-horizontal) aspect 
ratio. Experiments with fluid in a rotating annulus, using a forcing mecha- 
nism, enable to produce a zonal (azimuthal) jet subject to the Coriolis force. 
In such fast-rotating tanks equipped with ridges at the bottom (mimicking 
topography), flow patterns identified as 'zonal' and 'blocked' states are ob- 
served. In addition, transitions between the two states are found in a certain 
range of forcings (tank rotation and pumping rate) [15]. In the blocked state, 
streamlines tend to follow topography contours. 

This bistability is reminiscent of the phenomenon of atmospheric 'block- 
ing': On interannual time scales, large anticyclones form in the Northern 
Hemisphere, blocking and deflecting the nearly zonal flow (following lati- 
tude circles) [18]. Analogous configurations are observed in the north Pacific 
ocean, where the Kuroshio Extension forms an eastward mid-basin jet. The 
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Kuroshio is seen to oscillate between an intense jet-like (zonal) state and a 
weaker meandering (blocked) state. 

Our intuition is that the qualitatively different states predicted by our 
phase diagrams could be related to the different regimes observed in geo- 
physical flows, providing that realistic geometries are considered (annular 
domain for example, coastline geometry, bottom topography). Fluctuations 
would be responsible for the transitions between different equilibria. 
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5 Appendix: Grand canonical solutions 

Let M g = {q'\(q'ei) = 0}. We may decompose h = h\e\ + h', with h' G M g . 
Through Atp = q — h, we have ip = — ^^^ ei + tp', with ip' G M g . Let us 
denote the Gibbs free energy (4) with 7 = by Q . 

g [ q } = g [A,q'}= [ 

Jv 

The second-order variation of Qq with respect to q' is 

8 2 g [A,q'} = [ s"{Ae x + q')5f -(if S^'Sq'. (13) 

It is straightforward to prove a generalization of the Poincare inequality in 
the subspace M g (any q' G M g may be written q' = J2i>2 Qi e i)j which yields, 
for (3 < 0, the inequality 

5*g [A,q']>U + £j Jjq' 2 , 

where s 9 A := min re x>{min g / s"{Ae\ + q'{r))}. Therefore, if (3 > —s 9 A \2, Go is 
convex with respect to q 1 and we denote by q' eq the unique solution to the 
minimization problem 

G (A)=mmg [A,q f ]= [ 
1 Jv 

2 Ai 

For > —s 9 A \2, q' eq is the unique critical point of Q$ with respect to q' . 
It satisfies 

J {s\Ae 1 + q' eq ) - (3^ eq )5q' = for all 5q' G M g , 

therefore there exists a } 6l such that 

s'{A ei + q' eq ) - (3ip' eq = a g e x . (15) 

We compute the solution to (15) perturbatively around (A, q') = (0,0), in 
order to obtain an asymptotic expansion for Gq (A) around A = 0, and hence 
determine the type of phase transitions to expect in the vicinity of (f3 < 
—s" n \i , 04 = 0). Remark that if ci2n < for all n > 2, then s'^ = min g {l — 



s(Ae 1 + q')-^\q' 



ti) 



+ 



2 Ai 



s(Ae 1 + q' eq )-^' eq (q' eq -h') 



+ 
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3a 4 g 2 — o(g 3 )} = 1 and s'^ — > 1 as (a 4 , g) — >■ (0, 0); s(q) is said to be strongly 
convex 2 . 

We have the Taylor expansion s'(q) = q — a^q 3 — a^q 5 + o(g 6 ). Substi- 
tuting this expression into (15), and projecting (15) orthogonally onto M g 
(projection of x being denoted by P(x) := x — (xei)ei), we get 

q' eq - a 4 A 3 P(e 3 ) - 3a 4 A 2 P{e\q' eq ) - 3a 4 AP {e x q%) - a 4 P(q' 3 q ) + 
+ 0{A\ A% q , A 3 q% A 2 q' 3 q , Aq' e % q%) = ^ eq . 

At lowest order in the asymptotic expansion of q' eq (q' eq = q' + higher powers 
of A), we have 

% - Wo = A"0o ~ 0^0 = a 4 A 3 P(el); 
the linear operator : ij)' i— > Aip' — (3ip' is invertible in the subspace M g for 
/3 in the vicinity of —Ax- Thus, we get 

^ = a A A 3 C^P(e 3 ) =: ^A 3 , 
q' = a A A 3 AC/P(e 3 ) =: q' A 3 . 

Now, we compute the asymptotic expansion of Go (A) using this perturbative 
result: substituting q' = q' A 3 + o(A 3 ) into (14), we get 



Go(A) = -Z-hiA + 2 ( 1 + t M " "f < e I>^ + ( 17 ) 



A 6 + o(A 6 ). 



The parity of Go (A) is broken by h x ^ 0. Let us take hi = until fur- 
ther notice. Note that up to quartic order, only mode e x contributes — the 
perturbation q' contributes only from order 6 and up. 



6 Appendix: Canonical solutions 

Recall 

F{/3, r) = min ^F[q] = jf s{q) + (3£[q] \ r[q] = T 
We reduce the set of independent variables to {qi}i>2, as in [16]. Since 

r = X)ft< e *)> ( 18 ) 



2 http : //en . wikipedia . org/ w/ index . php?title=Convex_f unction&oldid= 
466661785#Strongly_convex_f unctions Accessed January 5, 2012. 
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we may decompose 



r / (eO \ r 

- T ei+> e^-— — ei =:— — e! + g c , 19 

(ei) V (ei> / (ei) 



so as to consider the minimization of T with respect to q c : 

r 

ei + q c 



min {J 7 ^] | r[q\ = V} — min < T 



(ei) 



(20) 



Let M c = {g c , determined by {Qi}i>2\{Qc) = 0}. M c is asubspace complemen- 
tary to the line spanned by e±. M c is a subspace orthogonal to 1. Canonical 
solutions live in M c . 

Note 

r 



Sr ■= mm < min s" I - — -e\ + q ( 

reV [g c SM c \ ( e l) 

and consider 5 2 J-", the second-order variation of J 7 with respect to q c . In 
appendix 7, we prove a generalization of the Poincare inequality in the 
subspace M c , leading to S 1 ? > (sr + (3/ X c )(Sq^) . We have introduced 
A c := minjA*, A^}, corresponding to the vanishing of the quadratic part of 
J-"[Tei/ (ei) + q c ] (denoted by Qjr) at (3 = — A c . The space M c is a direct sum 
of the subspace generated by eigenmodes of zero domain average, {e^}j>i, 
and the subspace generated by all the other modes. In the former subspace, 
Qjr vanishes at (3 = — X[ along e[. In the latter subspace, Qjr vanishes at 
(3 = — A* along e*, where A* is the smallest value of —(3 such that 



*) = -Ern = ' (21) 

i>l 



K + (3 



The interested reader can find details about the above function in [7]. Any- 
how, there are no phase transitions in the canonical ensemble for (3 > — SrX c - 

In appendix 8, we detail the computation of the T = solutions. For 
circulation T, the expression of the 'linear' solution (04 = 0) is 

Q (B > -A n = — V ^Mi e" 

M = -Ac, r) = -S- Y: §^- b 4 ± Ae c . (22) 

We can see that a nonzero circulation will introduce a symmetry breaking 
into the normal form (31)-(32). We consider a small circulation |T|, for the 



26 



description to remain close to the zero-circulation case. Also, this is required 
by the low-energy limit and the vicinity of [3 = (3 C . Because of the A y —A 
symmetry breaking, due to T ^ 0, the second-order phase transition vanishes, 
leaving a phase diagram with a critical point. 



7 Appendix: Poincare inequality in the canon- 
ical ensemble 

In this appendix, we prove a generalization of the Poincare inequality to the 
case with fixed circulation, i.e., in the subspace M c . 

Let q G M c . Then, 5q = Y.i>2 ^( e * _ M e 
and 5^ = - £ 4 > 2 5 Qi (g - . We have 

Jv i>2 \ ei) i,j>2 



i>2 x v x/ i,j>2 



Now, 



is positive definite if and only if /3 > — minjA^, A*} = — A c . 

Since -/3 j v 54>5q > -13/13 J v 5q 2 for all G [-A c ,0[, then the best 
(greatest) lower bound that we can obtain is 

-13 f 8^8q>(- f 5q 2 . 

JV A c JV 

8 Appendix: Lyapunov— Schmidt reduction in 
the canonical ensemble 

In this appendix, we derive the phase diagram for the canonical solutions at 
zero circulation. Consider the following canonical variational problem (let us 
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drop the D subscript in the integral notation): 

^ = {/ ^ ~ 2 / ^ ' I q = °} ' 
A critical point is q such that 

J s '(<?)<5<? - P J i>fiq = for all Sq E Q such that J 5q = 0, 
or, equivalently, using the Lagrange multiplier rule, 

f(q,T,P):=< f( R , f n (23) 

h(q,i\P) = / q = o, 



where 7 6 R is the Lagrange parameter associated with the conservation of 
(zero) circulation. 

The system (23) is to be solved in the variables (g, 7), while the bifurcation 
parameter is /3 £ R. Let us denote the variable by X = (q, 7) and the variable 
space by E. Please do not get this notation mixed up with the energy, which 
we never mention in this appendix. / maps £xR into E. For any £ M, 
we have the trivial solution X = 0. We want to determine the bifurcations, 
which the system may undergo, from this trivial solution. 

For a bifurcation to occur, the Jacobian matrix of (23) has to become 
singular, i.e., there must exist a nontrivial vector u c = (q c , j c ) £ E such that 
Dxf{^'i /3)[it c ] = for a certain j3 = j3 c . We have 

s"(0)5q c -(35^ c + 6 lc 
D x f(0;f3)[u c }= I f I eE, (24) 




with Aip c = q c . Let us endow E with the scalar product (-|-), defined as 
follows: for X k = (q k ,j k ) £ E, k = {1,2}, 

(Xi|X 2 ) = {q x q 2 ) + 7! 72 = J (qiq 2 ) + 7i72- (25) 

A complete orthonormal basis for E is {uj}j>o, where uq = (0, 1) and itj = 
(ej,0) for i > 1. Dxf{0;/3) is self-adjoint since 



(X^DxfiO; 0)[X 2 \) = j qi (s"(0)5q 2 - /36^ 2 + 5 l2 ) + 71 J Sq 2 

= s"(0)(q 1 8q 2 ) - /3{qi6if) 2 ) + (?i^7 2 ) + (71^2) 
= (D x mP)[X 1 ]\X 2 ), 
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so Dxf{0] P) may be diagonalized in {ui}j>o, and its eigenvalues are real. 
The equalities (qi5q2) = (^Qi), {qxty-i) — (<?2^i); an d so on ) come from 
the Euclidean-ness of E. Indeed, let k = {1,2} and q k = YliQk^i- The 
tangent vector 5q k = ^Sq^ei is along q k , so for all % > 1, 5q kii = a k q k)i . 
Now, a± = ci2 because q\ and q2, belonging to the same space, must be 
mapped onto their tangent space with the same coefficient. Decomposing 
the variables in the Laplacian eigenbasis {ej}j>i, 

5q = ^2 / 5q i e i ] Sip = - ^ ~^ e i I h = h^^i)^ ( 26 ) 

i>l i>l 1 i>\ 

it is readily seen that u c is either along (e^,0) at /? = — s"(0)X' i , or along 
(e*, 1) at = — s"(0)A*. Indeed, we identify 



-fry ^ n fOT a11 i > 1, 

V(0)A. t + /3 



and we require 



/x \ sr x X ^ 2 5l f ( 13 \ 



where the / function is (21). We have noted 



e " 



i> 



A 



-e, ; . 



Note again that SX belongs to the tangent space of E, but E is Euclidean, 
so {X E E I SX = au c , aeR} = {XeE\X = au c , a E E}. 

The first bifurcation, and hence, phase transition, to occur is found at 
P c = — s"(0)A c (considering a decreasing ft). Let u c = jV(e*, 1) in case i), 
u c = (e[,0) in case ii). J\f is the normalization factor (((e*) 2 } + l) -1 / 2 . 

Let us denote the operator Dxf{0] Pc) by J. J maps E into E. Let E c 
be the null space (kernel) of J. It is the subspace generated by u c ; it is 1- 
dimensional in E (it is a line). Let us show that the range of J is orthogonal 
to E c . This is the case if and only if (y|it c ) = for any Y in the range of J. 

i) Let us show that 

((s"(0)5q - PM + 8j)e*) + (Sq) = 0. 
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Let ip* be the vector such that Aip* = e* and ip* = on dT>. We have 
(0e*) = (g^*) (straightforward when decomposing in the Laplacian 
eigenbasis), so 

((s"(0)5q - PM + h)e*) + (6q) = ((s"(0)e, - + l)Sq) = 0. 

(27) 

Indeed, the last parenthesed term is the first component of Ju c (Ju c = 
0). 

ii) Let us show that 

((s"(0)Sq-/3M + ^)e[) = 0. 
We have C = -s"(0)A' 1 , so 



{(s"(0)5q - PM + 



s»(0)e[ + |ei) 6q} = 0. 



Therefore, the kernel of J is orthogonal to the range of J. We can then apply 
classical bifurcation theorems [8]. Let E\ be the orthogonal complementary 
subspace to E c in E (E\ is the range of J). There exist X(A, 0) 6 E\ {A e M) 
such that X(0,(3 C ) = and §f (0, /3 C ) = 0, so that we may decompose the 
variable X as follows: 

X = X(A, (3) = Au c + X(A, (3) = Au c + (q(A, P),>y(A, (3)) . (28) 

We will also use the notation ip for the vector in Q such that Aip = q. 
Besides, there exists a projector Q : E — > Ei, QX = X — (X\u c )u c such that 
Qf(X; /3) = Qf(Au c + X;/3)=0 for all A,/3eR. Now, 

f(Au c + X;(3) = Qf(Au c + X;(3) + (f(Au c + X; (3)\u c )u c 

so the bifurcation problem (23) is equivalent to (reduces to) the scalar prob- 
lem 

h(A,/3) := (f(Au c + X(A,f3);f3)\u c ) = 

(Lyapunov-Schmidt reduction). From the normalization of u c , we have 
f(Au c + X; (3) = h(A, (3)u c . Explicitly, this writes 

i) 

s'(ANe* + q)- P{ANip* + $) + AM + 7 = jVh(A, f3)e*, 

(29) 



J AJVe* + q=JVh(A,f3); 
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s'iAe', + q) + 0(£e£ - ^) + 7 = MA /3)e'i, 

Al (30) 



J Ae[ + q = 0. 



We may notice that for (A, X, h) solution, (—A, —X, —h) is also a solution, 
so that h and X are odd in A. Therefore and |^ are also odd in A, and 
so on. 

We know that F is even in A. We have F(A = 0) = 0, so the lowest order 
of F is quadratic. We determine the successive coefficients (of each power of 
A) in F from its successive derivatives w.r.t. A, evaluated at A = 0. Because 
(qe c ) = 0, we also have (^f e c ) = 0, and so on with all the derivatives with 
respect to the scalar A. All these properties lead to drastic simplifications in 
the computation of gil(^4 = 0) and ^(s(A = 0), leaving us with 

i) 



F(A) = ^(0) + f-J A> - j^fjy^ + o(A% (31) 

") 

F(A) = \ (V(0) + £) A* - ?±(eM* + o{A% (32) 

Bifurcation-wise, it is shown that 

dh d 2 h 
h(0,p c ) = 0, — (0,/3 c ) = 0, —((),&) = 0, 

but 

-(0,/3 c )^0; S gn |-(0,/3 c ) =-sgn(a 4 ). 



Therefore, the bifurcation will be determined (qualitatively) by the cubic 
nonlinearity of h in A (corresponding to the quartic nonlinearity of F in A, 
in the present paper). The sign of a 4 , i.e., the parameter for the nonlinearity 
in the q — ip relationship, determines the type of bifurcation at play: 

- If 04 < 0, the pitchfork bifurcation is supercritical, giving a second- 
order phase transition. 

- If 04 > 0, the pitchfork bifurcation is subcritical, giving a first-order 
phase transition (the higher-order nonlinearities yielding nontrivial branches 
beyond j3 = f3 c , at /3 < (3 C ). 
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